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Motivated by the very recent fabrication of sub-lO-nm-wide semiconducting graphene nanoribbons 
[X. Li et al, Science 319, 1229 (2008)], where some of their band gaps extracted from transport 
measurements were closely fitted to density functional theory predictions for magnetic ordering 
along zigzag edges that is responsible for the insulating ground state, we compute current-voltage 
(I-V) characteristics of finite-length zigzag graphene nanoribbons (ZGNR) attached to metallic 
contacts. The transport properties of such devices, at source-drain bias voltages beyond the linear 
response regime, are obtained using the nonequilibrium Green function formalism combined with the 
mean-field version of the Hubbard model fitted to reproduce the local spin density approximation 
description of magnetic ordering. Our results indicate that magnetic ordering and the corresponding 
band gap in ZGNR can be completely eliminated by passing large enough DC current through it. 
The threshold voltage for the onset of band gap collapse depends on the ZGNR length and the 
contact transparency. If the contact resistance is adjusted to experimentally measured value of 
~ 60 kf2, the threshold voltage for sub- 10-nm- wide ZGNR with inter-contact distance of ~ 7 nm 
is ~ 0.4 V. For some device setups, including 60 kSl contacts, the room temperature I-V curves 
demonstrate step-like current increase by an order of magnitude at the threshold voltage, and can 
exhibit a hysteretic behavior as well. On the other hand, poorly transmitting contacts can completely 
eliminate abrupt jump in the /- 1/ characteristics. The threshold voltage increases with the ZGNR 
length (e.g., reaching « 0.8 V for ~ 13 nm long ZGNR) which provides possible explanation of 
why the recent experiments [Wang et al., Phys. Rev. Lett. 100, 206803 (2008)] on ~ 100 nm 
long GNR field-effect transistors with bias voltage < 1 V did not detect the /- V curve signatures 
of the band gap collapse. Thus, observation of predicted abrupt jump in the I-V curve of two- 
terminal devices with short ZGNR channel and transparent metallic contacts will confirm its zigzag 
edge magnetic ordering via all- electrical measurements, as well as a current-fiow-driven magnetic- 
insulator-nonmagnetic-metal nonequilibrium phase transition. 

PACS numbers: 73.63.-b, 75. 75. -fa, 73.20.-r, 85.35.-p 



I. INTRODUCTION 

The recent surprising discovery of graphene^ — a one- 
atom-thick layer of graphite — has introduced in a short 
period of time a plethora of new concepts in condensed 
matter physics and nanotechnology, despite apparent 
simplicity of the two-dimensional honeycomb lattice of 
carbon atoms that underlies much of its unusual physics 
revolving around Dirac-like low-energy electronic exci- 
tations.- Examples include anomalous versions of meso- 
scopic transport effects,^ topological insulatorSf§ and 
low-dimensional carbon-based magnetism^i^ to name 
just a few. Since driving a system out of equilibrium 
typically corrupts its quantum coherence and suppresses 
quantum interference effects, basic research experiments 
have mostly been focused on the linear response regime^ 

At the same time, vigorous pursuit of carbon nano- 
electronicSfii^ envisioned around gated planar graphene 
structures that promise to overcome some of the diffi- 
cultiesi encountered by carbon nanotubes, has led to 
increasing number of experimentally demonstrated top- 
gated graphene field-effect transistor (FET) concepts. In 
these setups, micron-size graphene sheet o^'^°d^ or sub- 
10-nm-wide graphene nanoribbons^^ were employed to 
demonstrate room temperature graphcnc-FET operation 
with ON/OFF current ratios^^ up to 10^, high carrier 
mobility in sheet-based FETs,^° large critical current 



densities,^ and operating frequency reaching-'^"'^ ~ 26 GHz. 

These experiments pose a challenge for theoretical and 
computational modeling since they drive graphene nanos- 
tructures into far-from-equilibrium regime due to the fi- 
nite applied bias voltage. The task is more demanding 
than typical linear response-based analysi a^^d'^ of poten- 
tial graphene devices due to the need to compute self- 
consistently developed potential and charge redistribu- 
tion within the system in nonequilibrium current carry- 
ing state in order to keep the gauge invariance^^ of the 
/- V characteristics intact. 

In addition, the description of experimental devices 
often requires to include much greater microscopic de- 
tailsi^iii than captured by simplified effective models 
that resemble relativistic Dirac Hamiltonian for mass- 
less fermionsHiiSi or its parent single 7r-orbital nearest- 
neighbor tight-binding Hamiltonian.— These include 
atomic (such as the presence of hydrogen atoms which 
passivate edge carbon atomai^) and electronic struc- 
ture (probed by the bias voltage defined energy window 
around the Dirac point), self-consistent charge transfer 
effectsSS. that depend on the environment of an atom 
(tight-binding models are blind with respect to the charge 
of the system), and possibly more intricate manifesta- 
tions2ii2^ of electron-electron interactions in quasi-one- 
dimensional graphene nanostructures. 

For example, unlike the sheets of bulk graphene which 
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can be viewed as a zero-gap semiconductor jiSill mea- 
sured ratios of currents in ON and OFF states /qn /Iqff 
for room-temperature GNRFET— reveal that all of the 
sub-lO-nm-wide nanoribbons underlying the device were 
semiconducting. Furthermore, some of the energy gaps 
extracted from the operation of GNRFETs were closely 
fitted by the density functional theory (DFT) predictions 
for magnetic insulating ground state of ZGNR whose 
band gap is inversely proportional to GNR width. This 
is in contrast to non-interacting continuum Dirao^i or 
tight-binding^'' models of ZGNR which find only metal- 
lic nanoribbons with no energy gap around the Fermi 
level. Even larger band gaps, predictedSi by more com- 
plicated (non-self-consistent) many-body GW treatment 
of putatively enhanced electron-electron interactions in 
very narrow ZGNRs, were not seen in these measure- 
ments. 

A GNR is created by cutting a graphene sheet along 
two parallel lines. The recently developed nanofabrica- 
tion techniques for sub-lO-nm-wide GNR include direct 
STM tip drawing— and chemical derivation.— The 
GNRs produced by the latter technique were used for 
/-F curve measurements in Refs. IT^ and[23l Their crys- 
tallographic orientations were not identified. However, 
the fact that the number of sp^ bonds per unit length 
to be cut by chemical derivation in ZGNR is less than 
the number of bonds in armchair GNR (AGNR) suggests 
that chemical derivation is more likely to produce ZGNR 
rather than AGNR. 

Both AGNR and ZGNR are predicted to be semicon- 
ducting where the origin of their band gap is different. 
The band gap in AGNR is the consequence of quantum 
confinement and increased hopping integral between the 
TT-orbitals on the atoms around the armchair edge caused 
by slight changes in atomic bonding length. On the 
other hand, the band gap in ZGNR is due to staggered 
sublattice potential arising due to non-zero spin polar- 
ization around the zigzag edges. 

Although not confirmed by direct probing (such as via 
sophisticated spin-polarized STM techniques able to de- 
tect magnetic moment of individual atoms^l), the pos- 
sibility of peculiar carbon-based s-p magnetism (in con- 
trast to conventional magnetism originating from d or / 
electrons^) has been known since the early studies^^ of 
edge localized states due to special topology'^" of zigzag 
edges. These states have partially flat (within one-third 
of the ID Brillouin zone) subband, thereby generating 
large peak in the density of states at the Dirac point (i.e., 
the Fermi energy of undoped graphene). This makes it 
possible to easily satisfy the Stoner criterion^* for mag- 
netic ordering when (even tiny^^) Coulomb interaction 
is taken into account, which is the most likely way to 
resolve the instability brought about by the high density 
of states at the Fermi level in the nonmagnetic ZGNRs. 
Furthermore, the study of ZGNR magnetism has recently 
emerged as one of the major topics of theoretical research 
on graphene, reignited in part by the DFT calculations 
within the local spin density approximation (LSDA) that 



have described properties of such ordering from first prin- 
ciples^OS^ 

In this equilibrium picture, the ground state electronic 
configuration of both infinite^ '^^i^^ and finite-length^i 
ZGNR is characterized by ferromagnetic ordering of spins 
at each zigzag edge, antiparallel spin orientation between 
the two edges, and antiferromagnetic coupling between 
the two edges. Such compensated ferrimagnetic order- 
ing within ZGNR free of defects has zero total magnetic 
moment. Since opposite spin states occupy different tri- 
angular sublattices of the honeycomb lattice, the corre- 
sponding staggered potential induces^ the energy gap. 
The gap is inversely proportional to the width of the 
ribbons because the potential in the middle of the rib- 
bon decreases as the width increases (the band gap van- 
ishes within the room-temperature thermal energy win- 
dow when the width of GNR reaches ~ 80 nm) 

These findings have also motivated numerous propos- 
als for applications of ZGNR and graphene nanoislands 
with zigzag edges in spintronicS(^i^i^i^i2i despite the fact 
that no true long-range ordering in one-dimension is ex- 
pected at finite temperatures (for example, at room tem- 
perature the range of magnetic ordering along the edge 
is quantified by the spin correlation length estimated^^ 
to be ~ 1 nm). Moreover, virtually all known manifesta- 
tions of edge magnetic ordering in ZGNR have been pre- 
dicted within the framework of equilibrium theorie o^'^^i^^ 
or linear response transport calculations''^ which assume 
vanishingly small bias voltage. Except for the study of 
its modification, and ultimately destruction, in idealized 
infinite ZGNR due to the passage of finite ballistic cur- 
rent, very little is known on how such magnetism will 
manifest in the transport properties of realistic devices 
where finite-length ZGNR is attached to metallic con- 
tacts^i^ and biased by finite voltage applied between 
electrodes, as is the case of experiments on GNRFETs 
reported in Refs. [H and [H- 

Here we describe the fate of the band gap in two- 
terminal ZGNR devices where finite bias voltage brings 
them into a nonequilibrium steady transport state. Our 
results predict that passing a large enough current along 
ZGNR results in the destruction of spin-polarization 
around zigzag edges. This, in turn, causes the col- 
lapse of magnetic-ordering-induced band gap, and hence 
can lead to an abrupt step in the I-V characteristics 
of ZGNR. Nevertheless, this fundamentally nonequi- 
librium effect was not observed in recent experimen- 
tal measurementsi^ of the GNRFET I-V character- 
istics. Therefore, the second principal goal of our 
study is to provide explicit prerequisites for the ex- 
perimental observation of current-fiow-driven collapse 
of spin-polarized state in ZGNRs and the correspond- 
ing magnetic-insulator-nonmagnetic-metal nonequilib- 
rium phase transition. We note that phenomenologi- 
cally similar voltage-driven HOMO-LUMO gap collapse 
in a molecule attached to two electrodes was predicted— 
when these two levels (broadened by quasiparticle scat- 
tering) hit the bias window simultaneously, which to- 
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gether with our findings emphasize possibihty of highly 
intricate phenomena due to the complexity of nonequi- 
librium steady state in the finite-bias transport regime. 

The paper is organized as follows. SectionlTTlintroduces 
the Hubbard model in the mean-field approximation as 
a two-parameter fit to ab initio LSDA. In Sec. IIIII we 
describe the Newton-Raphson method used to accelerate 
the convergence of self-consistent spin-resolved electron 
density in the nonequilibrium state. The minimal basis 
set and the Hubbard model make the relatively expen- 
sive Newton-Raphson method much simpler and very ef- 
ficient when combining with the nonequilibrium Green 
function (NEGF) techniques. Section IIVI considers the 
influence of finite temperature and edge disorder on the 
spin-polarization of ZGNR in equilibrium. In Sec. |Vl we 
present our principal results: (i) the threshold voltage 
required to destroy the edge spin polarization increases 
with the ZGNR length, reaching « 0.4 V and « 0.8 V 
for ZGNRs of length ~ 7 nm and ~ 13 nm, respectively 
(Figs. [5] and [H]); (ii) since larger threshold voltages and 
higher turn-on current may destroy ZGNR, we propose 
that the length of ZGNR intended for the band gap col- 
lapse measurements and, therefore, /- V curve probing of 
the underlying magnetic ordering, should be of the order 
of ^ 10 nm. We also discuss in Sec.|V]the influence of the 
contact quality on the observability of predicted features 
of the /- V characteristics of ZGNR sandwiched between 
two metallic electrodes. We conclude in Sec. IVIl while 
providing technical details of the self-consistent electron 
density calculations in the nonequilibrium state in Ap- 
pendix 1^ 

II. ZGNR EFFECTIVE MINIMAL-BASIS-SET 
SELF-CONSISTENT HAMILTONIAN AS A 
TWO-PARAMETER FIT TO LSDA 

The texture of magnetic ordering in confined graphene 
nanostructures, with at least few carbon atoms'^* form- 
ing a zigzag edge, has been described quantitatively ei- 
ther by using the mean-field approximation of the Hub- 
bard (MFAH) model with single 7r-orbital per site^i^i^^ 
or DFT within different approximation schemes for its 
exchange-correlation density functional (such as LSDA)^ 
GGA,5 or hybrid BSLYPi^). Ahhough DFT goes beyond 
strictly on-site treatment of the electron-electron interac- 
tion U and the nearest-neighbor hopping t in the MFAH 
model, the parameters of the latter— 

i?MFAH = XI H i^a^i'^ + Cj<T^i'^) 

+ {4^iT ("ii - y) + ^uCij, (nit - y) } 

i 

+E E ^■4A'r, (1) 

i 'y=hl 

can be estimated by fittingly, the spin-unrestricted DFT 
band structure near the Fermi energy Ep = with that 



obtained from the MFAH Hamiltonian ([T]) defined on 
iVx-ZGNR honeycomb lattice. The values for t and U 
estimated in this fashion slightly depend on the choice 
of the exchange-correlation functional employed within 
DFT approximation schemes^i^ In Eq. ([T|), operator c[ 
(ci) creates (annihilates) an electron in the 7r-orbital lo- 
cated at site i = {ix,iy) of the honeycomb lattice. The 
third term, which is zero for charge neutral systems, ac- 
counts for the shift of the on-site energies due to Coulomb 
interaction with the applied electric fields or uncompen- 
sated charges in the system — the coefficients Vi are to be 
computed self-consistently in a standard DFT-like fash- 
ion, as elaborated in Sec. IIIII 

As customary)^ the width of iV^-ZGNR lattice is mea- 
sured using the number of zigzag longitudinal chains. 
The number of atoms comprising a single longitudi- 
nal zigzag measures its length. In the units of graphene 
lattice constant a = 2.46 A, the average width of ZGNR 
is I^ = aV^{N, - l)/2 and its length is L = a(A^^ - l)/2. 

The spin-resolved {a J, along the 2;-axis orthogonal 
to ZGNR plane) electron density on carbon atom at site 
i is given by the statistical expectation value 

riia = (clAa), (2) 

so that particle density at the same site is the sum 

^i^riit+n^. (3) 

These quantities have to be computed via the self- 
consistent loop,^° cither from the cigcnstates of equilib- 
rium systems^ or from NEGFs (Sec. IIII[) when finite- 
length ZGNR is attached to electrodes'^^ and brought 
into nonequilibrium state by the applied bias voltage. 
Once the self-consistency criterion is satisfied, the spa- 
tial distribution of magnetization density within ZGNR 
is obtained from 

gt^sSi = HB{n^ -nil), (4) 

where S? is the spin density and fj,B is the Bohr magne- 
ton. 

We chose to combine local orbital basis Hamilto- 
nian ^ with NEGF because it allows us to substantially 
accelerate self-consistent calculations in the nonequilib- 
rium current-carrying state (as discussed in Sec. lIIip . Al- 
though Eq. ID) is typically obtained through mean-field 
decoupling scheme^^- by starting from the full many-body 
Hubbard model for lattice fermions, it can also be justi- 
fied with the framework of LSDA. Furthermore, the latter 
provides simple and clear explanation of the expression 
for the total energy, which will be required for thermo- 
dynamic analysis of Sec. IIVI 

By using the spin-restricted self-consistent 
environment-dependent tight-binding model (SC- 
EDTB), which is specifically tailored to simulate eigen- 
value spectra, electron densities and Coulomb potential 
distributions for carbon-hydrogen systemS )^^'^^ we can 
establish the relationship between Hamiltonian ([T|) and 
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its LSDA counterpart^^ 



LSDA 



P" 

2m 



|r — r\ 
+ 14c([n];r) + Ay-, 



^ + Fpp(r) + V;,t(r) 



(5) 
(6) 



The first five terms in this one-electron Hamiltonian are: 
kinetic energy operator, classical Hartree potential, pseu- 
dopotential associated with core electrons, external po- 
tential, and spin-restricted part of exchange-correlation 
potential, respectively. They do not depend on the spin 
polarization and can be accounted by the SC-EDTB 
model. The EDTB aspect^ of the model assumes that 
hopping matrix elements of the tight-binding Hamilto- 
nian depend not only on the distance between the two 
atoms on which the basis functions are centered, but 
also on the arrangement of neighboring atoms (i.e., it is 
analogous to a DFT scheme that accounts for three- and 
four-centre integrals, and with atomic orbitals adjusted 
to atomic environment). The SC-EDTB model adds pa- 
rameters to this non-self-consistent EDTB part in order 
to describe hydrocarbon bonds while taking into account 
the self-consistent'^° charge transfer."'^';-^ The last term 
AV^^ in Eq. ([5]) is different for f and | spins, where the 
spin-dependent V^"^ ([n^, n^]; r) exchange-correlation po- 
tential in LSDA is^^ 



n„=n„(r) 



(7) 

The exchange-correlation energy per particle e^dni ,ni) 
is extracted^ from an electron gas with uniform densities 

Let us compute the on-site 7r-orbital matrix element 
for AV^'J. We borrow the orbital parameters from SC- 
EDTB, particularly 



(8) 



which is one of the four localized orbitals per carbon atom 
comprising the basis set. Here Up = 1.6085 A~^/^ and a 
is the normalization factor. By defining local relative 
spin polarization 



C(r) 



n^(r) — n^{r) 
7iT(r) + nHr)' 



(9) 



and by assuming that within the orbital range C does 
not depend on r, the spin-resolved electron densities for 
TT-orbital are 



(r) 



nF^' (r) 



:1-C 



(10) 



Equations ([7]) and ([8]) determine the matrix element 
{^%\AV^^Yl!l) as a function of C (*^ = *zXa where Xa 
is the spinor part of the wave function). 
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FIG. 1: (Color online) The spin-dependent contribution to 
the on-site matrix elements of LSDA and MFAH Hamilto- 
nians as the function of relative spin polarization C,. The 
dashed red line plots the expectation value (^'I|A\4L|^'J) for 
Perdew-Zunger parametrizatior4^ and SC-EDTB orbital pa- 
rameters, where electron density contribution from other sp^ 
orbitals to AT^L is neglected. The solid blue line plots spin- 
dependent contribution J^mfahIC] to the on-site matrix ele- 
ment of MFAH Hamiltonian with U = 2.7 eV. The two ver- 
tical lines indicate variation range for the spin polarization 
parameter C, within ZGNR. The right inset plots the loga- 
rithm of SC-EDTB sp^ DOS as a function of energy {Ef = 
is the Fermi level) . The left inset plots spin-dependent contri- 
bution E (dashed red) to the total energy Eq. (fT^ in LSDA 
and — (7 niiTiii (solid blue) contribution to the total energy 
Eq. (fTH)) in MFAH model as a function of 



Similarly, we can obtain the spin-dependent contribu- 
tion to on-site matrix elements of MFAH Hamiltonian ^ 



(11) 



which is a linear function of C,- As we demonstrate below, 
even in the strong electric field the relative change of the 
total TT-orbital electron population does not exceed 1%. 
This means that ni in Eq. (Ill|) can be assumed equal to 
unity. 

Figure [1] plots the spin-dependent part of the on- 
site matrix elements of LSDA and MFAH Hamiltonians. 
Since in ZGNR systems the polarization C, varies within 
the interval [—0.3, +0.3], we find a good fit between ma- 
trix elements of MFAH and LSDA Hamiltonians in this 
range, thereby justifying the usage of Eq. ^ instead of 
more complicated Eq. The fit also sets the param- 
eters t ^ U = 2.7 cV for the MFAH model. The right 
inset in Fig. [T] plots the SC-EDTB contribution of sp^ 
carbon orbitals to the density of states (DOS) for a typ- 
ical single-layer graphene system composed of GNRs of 
different types and widths. As follows from the inset, 
the orbitals other than do not have any contribution 
to the DOS within [—4 eV, -f 3 eV] interval around the 
Fermi energy. This suggests that usage of single 7r-orbital 
per honeycomb lattice site in Hamiltonian ([T]) should be 
sufficient in the cases when the applied bias voltage does 
not exceed ±2 V. 

From this analysis, as well as from the mappings^ii^ 
of DFT calculations to simpler MFAH model or the fact 
that DFT results obey the Lieb theorem.*"' for the exact 



5 



ground state of the Hubbard model on charge neutral 
bipartite lattices, we can conclude that second neighbor 
hopping and intersite Coulomb repulsion (present in the 
DFT calculations) do not modify the relation between 
lattice imbalance and total spin of the ground state war- 
ranted for the Hubbard model for which these couplings 
are absent. Thus, given that the replacement of LSDA by 
MFAH Hamiltonian is reasonably justified, we proceed to 
derive expression for the total energy as a function of 
based on the solution of Eq. ([1]). The total energy within 
the LSDA framework is 



Kj 



,3 n{r)n{r') 



2 



r — r' 



(12a) 



+ y d^rn(r)exc(nT(r),nx(r)). (12b) 

Here are the eigenvalues of Eq.®, /(e) is the Fermi 
function, n is the chemical potential [chosen to satisfy 
/ d^rnir) = N where N is the total number of elec- 
trons], Zi are atomic core charges, and are nuclear 
coordinates. The self-consistent computation of spin and 
particle densities shows that transition between the spin- 
polarized and the spin-restricted solution does not result 
in any noticeable change in the total electron density. We 
also assume the same atomic coordinates for the entire 
range of interest for C,. Therefore, the second and the 
third term in Eq. (jl2ap do not depend on spin polariza- 
tion. If we assume that C is uniform, i.e., independent 
of r within single orbital range, then by substituting the 
electron densities Eq. (fTU]) into E[n,C] we obtain nearly 
quadratic dependence of E on as shown in the left inset 
(dashed hue) of Fig. [TJ 

The expression for the total energy within the MFAH 
model is given by^ 



TT^Total r ^^ 



E ■^(^^ ~ '^^^^ ^ ^E "-iT^-U- 



(13) 



The second term in Eq. ^3]) is plotted (solid line) as 
a function of C in the left inset of Fig. [1] and represents 
approximation for E in Eq. p2p plotted in the same inset 
(dashed line). As demonstrated by Fig. [I] the LSDA 
matrix element averaged over the range C G [~0-3, 0.3] 
is positive due to slightly superlinear dependence on C,. 
At the same time the average of the Hubbard matrix 
element over the same interval is exactly zero due to its 
linear dependence on C^. That is, on average MFAH model 
underestimates the on-site Hamiltonian matrix elements, 
but it overestimates E. The partial error compensation 
in Eq. ([T^ makes U = 2.7 eV a reasonable choice for 



approximation of both the LSDA single particle energies 
and the total energy by a simpler MFAH model. 

The energy expression also provides a transpar- 
ent explanation for the origin of magnetic ordering in 
ZGNRs. The value of C = — 1 in Fig. [1] corresponds 
to spin-t electron surrounded by spin-| electron density, 
while C = 1 is associated with spin-f electron surrounded 
by spin-f electron density. Therefore, the Hamiltonian 
matrix elements favor the spin polarization. At the same 
time, the second term in Eq. (jl3p favors the non-spin- 
polarized solution (see left inset in Fig. [1]), so that the 
competition between this term proportional to and the 
band energy proportional to C determines the appearance 
of non-zero spin polarization. 



III. NEGF WITH ACCELERATED 
CONVERGENCE SELF-CONSISTENT 
ALGORITHM 

We employ the NEGF formalism^ for the computa- 
tion of nonequilibrium spin-resolved electron densities by 
starting from the MFAH Hamiltonian Eq. ([1]). Assum- 
ing a two-terminal system, composed of finite-size ZGNR 
attached via semi-infinite ideal leads to the left (i) and 
right (i?) macroscopic reservoirs where electrons thermal- 
ize to be characterized by the electrochemical potentials 
> the nonequilibrium electron density 



n = diag [D] 



(14) 



in the phase-coherent approximation (i.e., in the absence 
of dephasing and inelastic processes^) is obtained from 
the following density matrix 

D = -i J dElm[G{E)]f{E-fin) 



+ CXD 



dERe {G{E)lm['EL(E)]GHE)} 

x[f{E~^^L)~fiE~^,n)], (15) 

Here G is the retarded Green function matrix and 
diag [. . .] returns vector composed of the diagonal ele- 
ments of its argument. Because Eq. ([1]) assumes zero 
overlap between the orbitals, only diagonal elements of D 
contribute to electron density in Eq. (fH|) . The retarded 
self-energy matrix Sl is introduced by the "interaction" 
with the left lead — it determines escape rates of electrons 
into the left reservoir."*^ The density matrix in Eq. (jl5|) 
is split into equilibrium (first term) and nonequilibrium 
(second term) contributionsj^i^i^ taking into account 
that left-lead states are filled up to fi^ and right-lead 
states are filled up to energy /i^. Integration over energy 
is performed using the poles summation algorithm.^? 

For a small difference between /x^ and hr, the self- 
consistency can be achieved by applying the Broyden 
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convergence acceleration methodr^i^i^ which has two 
major advantages. First, the Broyden method is com- 
patible with the recursive algorithm for construction of 
the Green functions and self-energies, where recursion 
is extended to allow for the computation of local quan- 
tities inside the sample rather than usual transmission 
function and conductance determined by it; 53,54,55 -j-j^g 
simplest version of such algorithms starts by partitioning 
the quasi-linear system into slices (described by a much 
smaller Hamiltonian matrix) in a such way that only the 
coupling between the nearest neighbor slices is present. 
Then, the recursive algorithm is applied to propagate 
the self-energies from the contacts and to build the Green 
functions for each slice. The nonequilibrium electron den- 
sity for each slice is derived locally from the Green func- 
tions and the self-energies for the given slice. The com- 
putation time scales linearly with the number of slices 
and cubically with the size of the matrices (Hamiltonian, 
Green functions, and self-energies) associated with a sin- 
gle slice. Second, the Broyden method adds 0{N) ex- 
tra operations and hence does not slow down the single 
iteration. However, the reduction of the iteration num- 
ber achieved by the Broyden method is appreciable. For 
equilibrium problems considered in Sec. lIVI the number of 
iterations required to achieve 10^^" maximum difference 
between the input and output electron densities using 
the Broyden method is about 30, while the scalar charge 
mixing requires several hundreds of iterations. 

The Broyden method works well when the correlation 
between the electron density and the potential is local, 
i.e., when the local potential distortion results in a local 
self-consistent density change. Conversely, in the case of 
non-local correlations the Broyden method performance 
rapidly deteriorates. The nonequilibrium electron den- 
sity in the coherent ballistic approximation constitutes 
the perfect example when the Broyden method fails. The 
reason for this is that electron-potential correlations be- 
comes completely non-local: the change of the potential 
at one contact can shut off the electron flux through 
the entire system and cause the system-wide electron 
density redistribution. The "brute-force" alternative to 
the Broyden algorithm is the Newton-Raphson method, 
which slows down each iteration by an order of magni- 
tude, but reduces the number of iterations to less than a 
dozen and guarantees the convergence towards the self- 
consistent solution. 

Because the convergence under nonequilibrium con- 
ditions constitutes the major computational problem, 
we present the details of the NEGF-adapted Newton- 
Raphson method employed in our study. The first or- 
der Taylor expansion for the retarded Green function 
(S = Si + SflJ 



G{E) [£; + H + S] 



(16) 



with respect to the Hamiltonian variation (5H is 

SG{E) = [E - (H + sn) - ^]-'^ - [E - n - T,]-'^ 

= G-(5H-G. (17) 



At this point the minimal-basis-set of the MFAH model 
comes into play — according to Eq.([T|) only diagonal ma- 
trix elements are affected by the electron density distri- 
bution. In the following, 611 denotes the change of the 
Hamiltonian due to a small variation of the electron den- 
sity, which means that SH is a diagonal matrix. Vector 
6h denotes the diagonal elements of ^H. Using the sym- 
metry of the retarded Green function matrix associated 
with the real Hamiltonian, the variation of the electron 
density Sn with respect to Sh can be written as (all quan- 
tities depend on energy E which is omitted for brevity): 



Sn — A - Sh, 



(18a) 



+ 00 



A = 



J dEIm[GtE)G]fiE ^ hr) 

— OO 

- - / dERc {Gt ® (G • Impi] • Gt) } 

^ J —OO 

x[f{E-fiL)-f{E~fiR)]. (18b) 



Here the symbol (E" between matrices denotes element- 
wise product of two matrices, so that the element of, e.g., 
G (g) G is {GpqY'. The computational complexity of the 
integrand in Eq. (jl8p is O (-/Vq) per energy point, where 
Ng is the size of matrix G. 

In the spin-unrestricted case the electron density vec- 
tor is composed of n| and n| sub-vectors. Therefore, we 
can rewrite Eq. (jl8p as a matrix equation: 



5n\ 



T - 



AT 
A^ 



(19) 



Matrices A"^ are computed using Eq. (jl8bp with Green 
function matrices G'^'^ plugged in. The matrices G°'°' are 
obtained by inverting via Eq. pB]) the corresponding di- 
agonal block H'^ (for spin-cr electrons) of the matrix rep- 
resentation of Hamiltonian ([T|) . We assume that there is 
no spin polarization in the leads, so that the self-energies 
Si and TiR are the same for both spin polarizations. 

In the framework of MFAH model, the on-site poten- 
tial variation vector is the linear function of the density 
variation vector 



(5h| 
5^1 



Q Q + C/I 
Q + C/I Q 



(5n| 



(20) 



Here Q is the Coulomb interaction matrix computed for 
TT-orbital wave functions with SC-EDTB parameters in 
Eq. (HI) using standard DFT approach. The dot-product 
of the i*^ row of matrix Q and uncompensated 7r-orbital 
electron density (n — 1) plus the potential shift due to the 
external electric field equals the coefficient Vi in Eq. ([1]). 
The identity matrix I has the same dimensions as Q and 
H'^. By defining matrix B as 



B 



AT 
AT 



Q Q + f/I 
Q + C/I Q 



(21) 
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we can relate, in the first order approximation, the re- 
sponse of the output spin-resolved electron density with 
respect to a small variation of the input density: 

f^^) 

\ / out \ ^ / in 

If riin is the input density of the self-consistent loop, 
and Hout is the corresponding output density, the self- 
consistent solution can be written as 

nout + Suont = nout + B • Sn-m = ni„ + ^ni„. (23) 

Equation ((23l) allows us to compute Sriin for the next 
self-consistent iteration from riin and nout in the current 
iteration by solving the system of linear equations 

(Ib - B) • Suin = nout - ni„. (24) 

Here Ib is the identity matrix of the same dimension as 
matrix B. 

The main computational disadvantage of the Newton- 
Raphson method is that Eq. (jl8bp uses the full retarded 
Green function matrix G, rather than its diagonal part as 
does the Broyden method. This prohibits the usage of the 
recursive Green function algorithm, and requires to ap- 
ply the Newton- Raphson scheme to matrices containing 
the information about the entire system rather than to 
much smaller matrices containing the information about 
its slices. Given that the second term in Eq. (jl8bp must be 
evaluated for about one thousand different energy poles, 
we are limited to systems composed of relatively small 
number of carbon atoms — the largest out-of-equilibrium 
ZGNR-based two-terminal device treated in Sec. |V] con- 
tains about one thousand atoms. 



IV. EQUILIBRIUM THERMODYNAMICS OF 
ZGNR 

A. Finite-length ideal ZGNR 

With few exceptions, ^^'■^^ theoretical investigations of 
magnetic ordering in ZGNRs have concentrated largely 
on all-graphitic structures (with addition of different 
types of edge carbon atom passivation^^) . On the other 
hand, in experiments, the ultimate electronic contacts 
are metallic, as illustrated by sub- 10-nm- wide GNRFETs 
with Pd source and drain electrodesJ^i^^ We first an- 
alyze equilibrium magnetic properties of ZGNRs of fi- 
nite length, with no defects and bounded by perfectly 
formed zigzag edge, which are attached to metallic leads 
modeled as semi-infinite square lattice wires. The device 
setup is illustrated in Fig. (5] We assume that on the 
square tight-binding lattice of the leads only the nearest- 
neighbor hopping tsi = t ~ 2.7 eV is different from zero, 
and the coupling of the leads to central ZGNR sample is 
described by the same hopping parameter — t. 

At the Fermi energy [Ep = 0) of undoped graphene, 
such leads have maximum number of open transverse 




Legend; Q nt-n-l- - +0.3 Q nt-n^ - -0.3 



FIG. 2: (Color online) Bottom panel: Self-consistently 
computed equilibrium spin density within finite-length 
6-ZGNR attached to two semi-infinite square lattice leads at 
T = 293 K. Upper panel: The initial random spin density 
used to obtain the solution in the bottom panel. Thick red 
and thin blue circles denote spin-^ and spin-| densities, re- 
spectively, with the circle radius being proportional to spin 
density on the corresponding carbon atom. 



propagating modes, which can penetrate into ZGNR as 
evanescent modeSf2ii^ In fact, at clean armchair left and 
right interfaces of ZGNR mode mixing occurs, thereby 
effectively acting as disorder whose effect on lead- ZGNR 
contact transparency further depends on weather the 
lead is "lattice-matched" or "lattice-unmatched" to the 
honeycomb lattice of ZGNR.^^- The leads shown in Fig. [2] 
fall in the category of "lattice-unmatched" ones,"^^ as dis- 
cussed in more detail in Sec. IV CI 

The two-terminal device in Fig. [5] is macroscopically 
inhomogeneous, so that even in equilibrium it is more ef- 
ficient to use NEGF (with semi-infinite leads accounted 
through self-energies discussed in Sec. IIIip to obtain the 
texture of its spin polarization,'^^ rather than trying to 
match the eigenstates of the leads to the eigenstates of 
ZGNR. If one starts with the random spin polarization 
illustrated in the upper panel of Fig. [21 the self-consistent 
solution converges to the magnetization density forming a 
pattern of finite length segments whose spins are oriented 
in the same direction. For example, the lower panel of 
Fig. [2] displays one possible self-consistent solution orig- 
inating from the initial spin density in the upper panel. 
The magnetization texture within ZGNR lowers the total 
energy but at the same time it decreases the entropy by 
aligning electron spins. Thus, at finite temperature one 
can expect that spin density along zigzag edges would be 
represented by finite length segments of a given polariza- 
tion, as exemplified by the lower pane in Fig. [2] 

To estimate the length of a uniformly spin-polarized 
segment, we compute the free energy F for ZGNRs of 
finite length and find its minimum with respect to the 
number of such segments. This problem can be formu- 
lated as follows. Suppose there is continuous L unit- 
cells-long zigzag edge, which can have either uniform or 
fragmented spin polarization. The length Ig of the short- 
est possible fragment for any given value of Hubbard U is 
known from self-consistent calculations. Obviously, the 
edge cannot contain more than L/lg fragments. Every 
boundary between the two fragments adds additional en- 
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"ff Efi (meV) Q nT-n4. = -0.25 



FIG. 3: (Color online) Left panel: Free energy of the zigzag 
edge of length A''" = 101 (50 one- dimensional unit cells) as 
a function of the number of segments with uniform spin- 
polarization at two different temperatures. The minimum 
at Us « 2.2 indicates that at T = 273 K the ZGNR edge is 
partitioned into two segments, while the absence of the mini- 
mum at T = 250 K means that the number of partitions is less 
than one so that the whole edge is uniformly spin-polarized. 
Right panel: The average length of a uniformly spin-polarized 
segment as a function of the energy En associated with the 
"boundary" between two segments of opposite spin polariza- 
tion. 



ergy E-^i to the total energy of the edge, which is also 
determined from the self-consistent loop. If the number 
of segments with uniform spin polarization is < L/l^, 
there are "extra" L — nJs edge carbon atoms that can 
be distributed between Ug segments. Thus, this prob- 
lem maps onto a question: "In how many ways can we 
distribute L — rislg indistinguishable spheres among 
distinct boxes?" Its answer is simply 



[L - nJs + ns - 1)! 

{L - njsy.iris - ly: 



(25) 



By applying the Stirling approximation for large factori- 
als, Eq. (I^S)) can be transformed into 



^(L + n.(l-;.)-l)"=(^-'=)+^-^ 



(ns-l)2 ""{L-nslg) 



(26) 



so that the entropy related to the spin arrangement along 
the edges is S* = fcs In [W (n^, ls,L)]. The free energy is 
then given by 

F{n,,h,T,En) = 

(n, - 1) - keT In [W (n^Js, L)] . (27) 

The number of segments Ug in ZGNR which is L unit- 
cells-long at equilibrium is obtained from the condition 



OF 

dn. 



^ 0. 



(28) 



We use the self-consistent calculations, similar to the 
one displayed in Fig. [5] but for longer ZGNRs, to extract 



FIG. 4: (Color online) Equilibrium spin density within finite- 
length 8-ZGNR, whose upper zigzag edge is eroded, at 
T = 293 K. The ZGNR is attached to two semi-infinite square 
lattice leads of the same type as in Fig. [5] The dangling bonds 
of edge atoms are assumed to be passivated with hydrogen 
atoms (not show explicitly). The circle radius is proportional 
to spin polarization at a given atomic site. Thick red cir- 
cles mark spin-| density and thin blue circles are for spin-| 
density. 



the values for Ig and E-^. The average number of atoms 
substantially affected by the transition between the two 
zigzag edge segments of opposite spin polarization is ~ 4, 
as illustrated by the lower panel in Fig. O The energy 
of the "boundary" between two oppositely spin-polarized 
edge segments is computed from 



E 



E' 



Total 



n 



-^Total 



(29) 



The value of E-^i slightly depends on the ZGNR width 
and spin ordering type at low temperatures, which can 
beJ^ (i) antiferromagnetic (AF), when spin moments on 
one edge are antialigned to the spin moments on the op- 
posite edge; and (ii) ferromagnetic (FM), when spin mo- 
ments on carbon atoms on both edges point in the same 
direction. Both the AF and FM configurations of mag- 
netic moments have total energy lower than the nonmag- 
netic state. Moreover, the AF configuration is the ground 
state in narrow ribbons, while the energy difference be- 
tween the AF and FM states diminishes with increasing 
ZGNR widths Here i^Totai is the total energy of ZGNR 
with arbitrary fragmented sections of uniform spin polar- 
ization along the zigzag edge, i?xotai total energy 
of the same ZGNR in the AF ground state, and Nj, is the 
number of transitions between "f and J, polarizations on 
both edges of ZGNR. 

The left panel of Fig. [3] plots F vs. the number of 
uniformly spin-polarized segments for ZGNR of length 
— 101 assuming Ig = i and E^ = 80 meV (i?|x can 
be estimated from the energy gap of ZGNR in magnetic 
insulating state, see Fig. (S]). The minimum for the free 
energy at T = 273 K indicates that the edges of ZGNR 
will be most likely partitioned into two uniformly spin- 
polarized segments but with antialigned spins between 
the two segments. On the other hand, at T = 250 K the 
free energy does not have minimum, meaning that the 
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whole zigzag edge is now uniformly spin-polarized and 
the corresponding ZGNR is in the AF configuration. In 
the limiting case of an infinitely long ZGNR {L ^ oo, 
rig — > c», and L/us — Zavr), Eq- (|28p simplifies to 



exp 



-^Ti \ _ ('avr + 1) 



kT 



1,-1 



(30) 



where Zavr is the average number of carbon edge atoms 
in the segment with uniform spin-polarization. The right 
panel in Fig. [3] plots l^vr as a function of energy E'n. For 



E^l =80 meV, this length is la 



28 at T = 273 K, 



while it increases by a factor of two at T = 250 K. 



B. Finite-length ZGNR with disordered edges 

Since edge imperfections are expected to disrupt'^^ the 
magnetic ordering within ZGNRs, we plot in Fig. |4] the 
self-consistent solution of NEGF-MFAH equations for 
ZGNR with vacancies along one of its zigzag edges. The 
removal of a single zigzag chain fragment from ZGNR 
edge results in almost complete loss of correlations be- 
tween magnetic moments of edge atoms belonging to dif- 
ferent chains. This means that the length of a segment of 
edge carbon atoms carrying magnetic moments aligned 
in the same direction is determined by either topological 
disorder or thermodynamic disorder, whichever has the 
smaller characteristic length. 



ZGNR IN NONEQUILIBRIUM 
STEADY-STATE 



A. Infinite ideal ZGNR 

As it has been previously suggested,'^^^ passing a suffi- 
ciently large current along infinitely-long translationally 
invariant ZGNR can destroy completely its edge mag- 
netic ordering. To estimate the source-drain bias voltage 



eVds 



(31) 



necessary to wash out the spin density, two separate 
Fermi levels fiL and (for the left- and right-moving 
electrons, respectively) have to be used. Since this sys- 
tem is infinite and homogeneous, the nonequilibrium 
spin-resolved electron density on carbon atom at site i 
can be computed simply by using its propagating Bloch 
modes 



m 

X [/(em,fc - 



ML 



) + fiem.k-m)], (32) 



where C;™(fc) is the value of the Bloch amplitude obtained 
by solving Eq. ([1]) and the sum over m goes over all bands. 



0) 



0.04 



0.00 



(a) 




0.0 0.2 


0.4 


0.2 0.4 






\V ) 






2 




(b) 
















"44 










-2 


0.42 eV / 




0.37 eV- 


1 


2 




(c) 




(f) 


S 
Ul 


















-2 


0.43 eV / 




0.3SeY/ 


1 



1.0 O.S 1 

k 



O.S 1.0 

k 



FIG. 5; (Color online) Total nonequilibrium energy [(a) and 
(d)] and band structure [(b), (c), (e), and (f)] of infinitely long 
ideal 6-ZGNR [(a)-(c)] and 32-ZGNR [(d)-(f)] as the function 
of applied bias voltage fiL — f^R = eVds- Panel (a) plots the 
total energy per unit cell of 6-ZGNR for magnetically ordered 
AF configuration (solid red) and nonmagnetic state (dashed 
blue) . The threshold voltage at which spin-polarized solution 
becomes unstable is labeled by Vj. For Vds > Vt the only 
self-consistent solution available is the non-polarized metallic 
state. The threshold voltage decreases for wider 32-ZGNR in 
panel (d). Panels (b), (c), (e), and (f) plot the band structure 
for the bias voltage slightly below (Vi"") and slightly above 
i^t'") the threshold voltage. Only two subbands in the vicin- 
ity of the Fermi level are shown for clarity, while other sub- 
bands experience only minor changes when the transition be- 
tween spin-polarized and non-polarized states of ZGNR talies 
place. 



The values of /zl and are determined from the charge 
neutrality condition 



(nil +nii) = Njjc, 



(33) 



where A^uc is the total number of electrons in the unit 
cell of ZGNR. 

Figure \5\ depicts the dependence of the total energy 
and band structure of 6-ZGNR and 32-ZGNR with re- 
spect to the difference between the Fermi levels of the 
left- and right-moving electrons. Panels (b)-(f) pertain 
to the nonequilibrium case, but are similar to equilibrium 
band structure and can be used to demonstrate the dif- 
ference between magnetically ordered AF configuration 
and nonmagnetic state of ZGNRs. The spin polariza- 
tion lowers the Hamiltonian eigenenergies, which results 
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in the down-shift of the highest filled band in panels (b) 
and (e) with respect to non-polarized band structures 
plotted in panels (c) and (f). Note that both in equi- 
librium^^ and nonequilibrium studied here, magnetic ef- 
fects on higher subbands are negligible, so that Fig. [5] 
shows only the two subbands around the Fermi energy 
corresponding to spin-polarized [Figs. Eljb) and [U^e)] or 
non-polarized [Figs.[5Jc) and EJf)] edge states. 

When the separation between /i^ and fin exceeds the 
certain threshold value, the abrupt change from spin- 
polarized AF configuration to non-polarized state occurs, 
as demonstrated by Figs. EJa) and EJd). The step-like 
transition can be explained as follows. Suppose the strip 
is in the AF spin-polarized state and fiL > fJ^R- As the 
bias voltage is increased, the left-moving electrons with 
majority spin orientation lying above fiji are depopu- 
lated (i.e., exit through the left contact without being 
replaced by electrons from the right contact), and the 
right-traveling states with minority spin orientation be- 
low /iL are populated. This reduces the spin density and, 
hence, decreases the energy gap in the subband structure 
plotted in Figs. [Hb) and ^e). When the device has 
reached a steady state, the result of these processes can 
be viewed as an effective repopulation of the electronic 
states in ZGNR, where electrons are "excited" from the 
valence band into the conduction band, with each such 
"excitation" depopulating a state in the valence band and 
populating a corresponding state with opposite spin in 
the conduction band. The reduced spin polarization de- 
creases the band gap, thereby facilitating "excitations" 
that depopulate more left-traveling states with major- 
ity spin orientation lying above j-iR and populate more 
right-traveling states with minority spin orientation ly- 
ing below fiL. When separation between /i^ and fj,R ex- 
ceeds threshold of « 0.4 eV, this feedback mechanism be- 
comes positive^ and the spin-polarized insulating state 
collapses to non-polarized metallic solution whose sub- 
band structure is shown in Figs.[5l^c) and EJf). 



B. Finite-length ZGNR attached to 
multiple-linear-chain leads 

To simulate metallic contacts on the top of ZGNR, 
similar to Pd contacts of experimental devices in Ref . [H, 
we assume that every carbon atom in the contact re- 
gions (shaded with light yellow color in Fig. [7]) is con- 
nected to a linear tight-binding chain with hopping pa- 
rameter between the chain atoms tic — 2.7 eV. Thus, the 
metallic electrodes of a two-terminal device are simulated 
with a large number of non-interacting semi-infinite lin- 
ear chains. The choice for such model is stipulated by the 
metallic character of linear chains and convenient way to 
simulate a top contact connected to a large number of 
atoms within ZGNR. Two different hopping parameters 
between the chains and the carbon atoms are employed: 
tc — 2.7 eV simulates highly transparent contact, while 
tc = 0.27 eV is chosen to correspond to the contact re- 




0.0 0.2 0.4 0.6 0.0 0.2 0.4 0.6 
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Energy (eV) Energy <.eV) 

FIG. 6: (Color online) The /- F characteristics [(a) and (d)] 
and transmission function [(b), (c), (e), and (f)] for 11-ZGNR 
(width « 2.1 nm; length « 6.6 nm) two-terminal device de- 
picted in Fig. [7] Left column panels are for perfect coupling 
tc — 2.7 eV between ZGNR and multiple-linear-chain leads, 
while right column panels use tc — 0.27 eV which sets the con- 
tact resistance to ~ 60 kf2 (as in experiments of Ref. . In 
panels (a) and (d), dashed red line denotes 7-1/ curves in the 
AF state of ZGNR, while solid blue line denotes I- V curves 
for the same ZGNR after its magnetic ordering is destroyed 
by the applied bias voltage. Solid red curve in panels (b) and 
(e) is transmission function at point B in panel (a) or point 
E in panel (d), respectively, for nonequilibrium AF spin con- 
figuration illustrated in Fig. [T] Solid blue curve in panels 
(c) and (f) is transmission function at point C in panel (a) 
or point F in panel (d), respectively, for nonmagnetic metal- 
lic state. In panels (b), (c), (e), and (f), thin solid black 
line illustrates energy window f{E — eVds/2) — f{E + eVds/'2) 
over which T{E, Vds) is integrated to get the current at corre- 
sponding points B, C, E, and F in panels (a) and (d), while 
thin dashed green line and thin solid red line are transmis- 
sion functions of infinitely long ideal 11-ZGNR at zero bias 
Vds ~ and threshold bias Vds ~ V^t"" (defined in the same 
way as in Fig. respectively. 



sistance of ~ 60 kfl measured for GNRFET devices in 
Ref. [H Also, the width ~ 2 nm of all ZGNR we exam- 
ine below is selected to fall in the range of experimentally 
fabricated sub-lO-nm-wide GNRsfiS, all of which have ex- 
hibited semiconducting behavior in the /- V characteris- 
tics measurements. Note that all-semiconducting nature 
of ultra narrow GNRs appears to be a key advantage over 
single- wall carbon nanotubes (which in the similar diam- 
eter range are typically a mixture of semiconducting and 
metallic onesSS.), as candidates for the envisaged carbon 
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FIG. 7: (Color online) Spatial profile of spin [panel (a)] and 
charge [panels (b) and (c)] density within 11-ZGNR (width 
~ 2.1 nm; length « 6.6 nm) two-terminal device at the thresh- 
old voltage Vds = 0.44 V and for perfect coupling {tc — 2.7 eV) 
to multiple-linear-chain leads attached in yellow-shaded con- 
tact regions. In panel (a), the maximum difference between 
spin-1" and spin- J, electron density is nij — nij ~ 0.24. The 
electron density profile in panel (b) corresponds to spatial 
spin distribution in panel (a) and point B in the /- V curve in 
Fig. [6ja) . Panel (c) shows spatial profile of electron density 
in the nonequilibrium state of ZGNR marked by point C in 
Fig. [6ja), whose spin-polarization is completely washed out. 

nanoelectronics .— 

For devices described by an effective single-particle 
Hamiltonian, such as Eq. ([1]), with mean-field treatment 
of interactions and no dephasing processes the current 
at finite bias voltage can be computed from the NEGF- 
based formula^^i^>ia 

4-00 

HVds) = ^ J dET{E,Vds)[f{E-fiL)~f{E~fiH)], 

— oo 

(34) 

which integrates the self-consistent transmission function 

Tr {Tr{E + eVds/2)GTL{E - eVds/2)G^} , (35) 

for electrons injected at energy E to propagate from the 
left to the right electrode under the source-drain applied 
bias voltage hl — fJ-R = eVds ■ The energy window for the 
integral in Eq. (jM)) is defined by the difference of Fermi 
functions f {E — fi^) — f (E — ^n) of macroscopic reservoirs 
into which semi-infinite leads terminate. This "window" , 
at room temperature T = 293 K and for selected bias 
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FIG. 8: (Color online) Electrostatic potential profile along the 
line drawn in the middle of ZGNR in Fig.[7{a) for: (a) applied 
bias voltage V^s = 0.44 V and perfect coupling {tc = 2.7 eV) 
of the contact regions to multiple-linear-chain leads, where 
dashed red line is for magnetically ordered AF state corre- 
sponding to point B in Fig. |6ja) and solid blue line corre- 
sponds to nonmagnetic metallic state marked by point C in 
Fig.[6ja); (b) coupling tc = 0.27 eV and the applied bias volt- 
age Vda = 0.44 V (dashed red) corresponding to spin-polarized 
state marked by point E in Fig. [6jd), or Vds = 0.5 V (solid 
blue) corresponding to nonmagnetic metallic state marked by 
point F in Fig. |6jd) . 

voltage Vds, is shown explicitly (thin black solid line) 
in Figs. El El and [13] for three different types of ZGNR 
devices, where it encloses the portion of T{E^ Vds) vs. E 
curve which is integrated to get the current /(Vds)- The 
matrix 

TlAE) = I {^lAE) - ^lA^)) (36) 

accounts for the level broadening due to the coupling to 
the leads. 

In addition to self-consistent computation of spin- 
resolved electron density, which is required both in equi- 
librium and nonequilibrium, evaluation of Eq. (|34p re- 
quires to compute also the self-consistently developed 
electric potential profilers- due to the passage of current. 
The profile enters into MFAH Hamiltonian ([1]) through 
Vi term. This ensures gauge invariance of /- V charac- 
teristics, i.e., its invariance with respect to the shift of 
electric potential everywhere by a constant The tech- 
nical issues in converging nonequilibrium charge densities 
through self-consistent loop are discussed in Appendix [Xl 

The /- V characteristics of two-terminal ZGNR devices 
with both transparent and « 60 kfi resistive contacts are 
shown in Fig. [51 In both devices, at around bias voltage 
Vf w 0.44 V, current jumps abruptly by an order of mag- 
nitude. To understand the origin of the jump, we plot 
the transmission function T{E^ Vds) in Fig. [Sljust before 
[panels (b) and (e)] and just after [panels (c) and (f )] the 
discontinuity has occurred. These plots reveal insulating 
state on the low voltage side Vds < where transmis- 
sion probability is exponentially suppressed within the 
gap region. On the other hand, finite transmission prob- 
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FIG. 9; (Color online) The (a) /- V characteristics, (b)-(c) 
transmission function, and (d) electrostatic potential profile 
for 11-ZGNR with the distance between the contacts twice as 
large (« 13.2 nm) as in the devices in Figs. [5] and [7] The 
coupling between carbon atoms in the contact region and 
multiple-linear-chain leads is perfect {tc = 2.7 eV). In panels 
(a) and (d), dashed red line is for AF state of ZGNR, while 
solid blue line denotes /- V curve (a) or potential profile (d) 
for ZGNR after its magnetic ordering is destroyed by the ap- 
plied bias voltage. Point B marks Vds = 0.09 V, while points 
C and D correspond to Vds ~ 0.77 V. The transmission func- 
tion (thick red) in panels (b) and (c) is computed at points 
B and C, respectively, where thin solid black line illustrates 
energy window f{E — eVds/'i) — f{E + eKis/2) over which 
T[E,Vds) is integrated to get the current at corresponding 
points B and C. Thin green line is the transmission function 
of an infinite ideal 11-ZGNR at zero applied voltage. Panel (d) 
plots the electrostatic potential profile along extended device 
ZGNR -I- contact-regions for spin-polarized AF configuration 
(dashed red) at Vds = 0.77 V (point C) and non-polarized 
metallic state (solid blue) at Vds ~ 0.77 V (point D). 



ability appears in the metallic state on the high voh- 
age side Vds > Vt of such voltage-driven nonequihbrium 
phase transition.^^ When the transparency of the con- 
tacts is reduced in Figs.[6Ue) and [HJf), the transmission 
function acquires sharp peaks due to quantum interfer- 
ence effects in the absence of dephasing (e.g., the ampli- 
tude of the resonant mode builds up when electron waves 
leaking from the quasi-bound state in the ZGNR channel 
cancel the incident waves and enhance the transmitted 
ones), which is akin to resonant transmission through 
double barrier structures4^ 

The spin density corresponding to point B in Fig. [S] is 
plotted in Fig. [TJa), demonstrating that insulating state 
for Vds < Vt is magnetically ordered in a similar fash- 
ion as in equilibrium. The AF configuration in nonequi- 
hbrium shows that a small amount of spin polarization 
is present in the middle of the ribbon, as is the case 



of equilibrium AF configuration whose edge states pene- 
trate deeper (when compared to FM configuration) into 
the bulk^S 

Further microscopic insight about the charge dynamics 
of ZGNR driven by finite bias voltage and current flow 
is revealed by the profiles of electron density in Fig. [7{b) 
for the insulating state and in Fig. [7]^c) for the metallic 
state, as well as by the electric potential profile for these 
two states plotted in Fig. [51 For example, in the insu- 
lating state for both transparent and resistive contacts 
ZGNR device, the potential profile in Fig. [5]^ a) and [Sfb) 
is linear, as expected for tunneling. When the band gap 
collapses the potential profile shows well-defined voltage 
drops near the contact regions and almost constant be- 
havior within ZGNR channel, as expected for ballistic 
conductor. The electric potential profiles along ZGNR in 
Fig. [8{a) can be directly related to spatial distribution of 
charges in Figs.[71[b) and[7]Jc). That is, increased charge 
density around the contact regions in Fig. [TJc) for non- 
magnetic metallic ZGNR is responsible for the voltage 
drop being confined (solid blue line) mostly around the 
contacts in Fig. [5Ja). 

Since abrupt current jump at the threshold voltage, as 
the most distinctive signature of voltage-driven nonequi- 
hbrium phase transition between insulating and metal- 
lic states of ZGNR, was not observed in recent exper- 
iments^^'^"^ on sub-lO-nm-wide ZGNR nanoribbons, we 
also investigate how the threshold voltage is affected as 
the length of the ZGNR channel increases. By dou- 
bling the inter-contact distance, from 6.6 nm in Fig.[S]to 
13.2 nm in Fig. [51 we find that threshold voltage increases 
from « 0.44 V at point B in Fig. [H to at least « 0.77 V 
at point C in Fig. [51 The current jump by an order of 
magnitude and electric potential profile in this case are 
still similar to ZGNR two-terminal device in Figs. [H^a) 
and [H due to the fact that metallic nonmagnetic ZGNR 
above the threshold voltage is attached to highly trans- 
parent contacts. Nonetheless, the hysteretic behavior at 
discontinuity in Fig. [Sl^a) , we predict for short ZGNR at- 
tached to electrodes via transparent contacts, is almost 
completely removed in longer devices. 

The decrease of the ZGNR band gap size is not a grad- 
ual process, but it is triggered when the bias voltage 
reaches a specific value. This is particularly transpar- 
ent in infinite ideal ZGNR (with no voltage drop along 
ZGNR) of Sec. IV Al where the reduction of the band gap 
and diminishing of edge magnetization density is initi- 
ated when the bias voltage window becomes equal to the 
band gap value. With further increase of the bias volt- 
age, the gap and spin density decay quickly to zero.— 
On the other hand, in finite-length ZGNR, the region of 
negligible transmission function T{E, Vds) in Fig. [51 in- 
creases from panel (b) to panel (c) with increasing Vds 
from point B to threshold voltage at around point C. 
This is due to electrostatic potential tilting of the local 
band structure, so that band gaps in different regions of 
the ZGNR cover different energy ranges inside the win- 
dow where we observe T(i?, Vds) ^ in Figs. Wis) and 
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Conduction Band Bottom 




FIG. 10: (Color online) Schematic explanation of how ZGNR 
length affects the threshold voltage for nonequilibrium phase 
transition between its magnetically ordered insulating and 
nonmagnetic metallic states. Top left inset depicts the 
graphene pattern cut out of a single sheet, which could be 
used to observe band gap collapse and destruction of mag- 
netic ordering in ZGNR by current flow, (a) Schematic plot 
of the local band gap along the dashed line in the top right 
inset at zero bias voltage Vds = 0. (b) Local band gap in 
the spin-polarized state for short ZGNR under nonequilib- 
rium conditions Vds 7^ 0. (c) Electrostatic potential proflle 
after the band gap collapse, (d) Local band gap in the spin- 
polarized state for long ZGNR under nonequilibrium condi- 
tions Vds 7^ 0. Shading is panels (a), (b), and (d) denotes the 
occupancy of electron states — the lighter the shade, the less 
is the occupation probability, where dark color corresponds 
to occupation probability one. 



(d). 

Although one can expect that the values of the volt- 
age at which nonequilibrium phase transition takes place 
will increases with the thickness of the tunnel barrier 
region introduced by magnetically ordered ZGNR, the 
values of the threshold voltage (or, more appropriately, 
a window of voltages taking into account hysteretic be- 
havior) one can expect for realistic two-terminal devices 
are non-trivial. For example, in abstract infinite ZGNR 
of Sec. IV^ this value is hmited^^ to « 0.4 V. On the 
other hand, in realistic two-terminal devices we find 
that the threshold voltage increase with increasing length 
of ZGNR, whose schematic explanation is provided by 




'q 



I 



i 



FIG. 11: (Color online) Color-coded amplitudes of the trans- 
mission matrix elements \tpq{E)\ connecting transverse prop- 
agating modes at Ep ~ lO^^t in the left and right square 
lattice leads which are "lattice-unmatched" (left column) or 
"lattice-matched" (right column) to 8-ZGNR with collapsed 
band gap. The ZGNR width = 8 and length A^^" = 37 
are the same as in the two-terminal device studied in Figs. 1121 
and [131 



Fig. [To] depicting the band energy diagrams for the nar- 
row ZGNR bridging the two contacts. The spin-polarized 
state becomes unstable when the occupancy of the elec- 
tron levels in the valence band decreases below, and the 
conduction band occupancy increases beyond the thresh- 
old level. Under nonequilibrium conditions, the change 
in occupancy in the short strips becomes possible due to 
the tunneling through the band gap [Fig. riOlfb)]. which 
results in the subsequent band gap collapse [Fig. [TOTc)]. 
In contrast, the tunneling rate through the band gap in 
long strips may be too low, and the required change in 
electron population cannot be achieved for a given Yds, 
as illustrated by Fig. [TOTd) . 



Finite-length ZGNR attached to square lattice 
leads 



Since Sec. IV Bl suggest that the magnitude of abrupt 
current jump at the threshold voltage will depend on the 
quality of contacts through which the finite-length ZGNR 
is attached to external circuit, in this section we exam- 
ine two-terminal devices whose metallic leads, modeled 
as semi-infinite square lattice wire, are attached laterally 
(rather than vertically on the top of ZGNR contact re- 
gion as in Sec. IVB)) . Several different ways of attaching 
square lattice leads to the honeycomb lattice of ZGNR 
have been explored in a variety of recent quantum trans- 
port studiesi^i^i^ For example, one can attach "lattice- 
matched"— or "lattice-unmatched"—!^ leads illustrated 
in Fig. [11] The former case is matched in the sense that 
the lattice constant of such lead is equal to carbon-carbon 
distance in graphenci^ 
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nl4-n4' - 4-0.01 
nTfn4- - -0.01 



FIG. 12: (Color online) Spatial profile of spin [panel (a)] 
and charge [panels (b) and (c)] density within finite-length 
8-ZGNR (width ~ 1.5 nm; length « 4.4 nm) connected to 
square lattice leads at the threshold voltage Vds ~ 0.4 V: 

(a) spin density in the state marked by point B in Fig. 1131 

(b) charge distribution in the same state as in panel (a); and 

(c) charge distribution in the non-polarized state marked by 
point C in Fig. [T3] 



When the hopping parameters t^i = tc — t are selected 
to be the same in the square lattice region (tgi), across the 
interface (tc), and in graphene {t = 2.7 eV), the trans- 
port across the interface is nominally ballistic, rather 
than through a tunnel barrier generated by reduced tc 
or mismatched tsi and t. Nevertheless, the conductance 
of lead-ZGNR-lead device can be greatly reduced if many 
propagating modes from metallic leads do not couple well 
to evanescent modes in GNR, or to evanescent modes plus 
a single propagating mode at the Fermi energy of ZGNR 
with collapsed band gap ^^i^° (which are the only avail- 
able modes to carry transport in narrow ribbons with 
large gap between Ep — and the second subband^°). 
Although evanescent modes are effectively enabling dop- 
ing of GNR by metallic contacts, which is enhanced in 
short and wide GNRsf^*^ the armchair transverse inter- 
face of ZGNR coupled to square lattice generates concur- 
rently mixing of transverse propagating modes.' '^■^'^^ This 
is illustrated by the non-zero off-diagonal elements of the 
transmission matrix 



t{E) = ^Tr{E) ■ G{E) ■ ^Tl{E) 



(37) 



in Fig. [TT] for metallic ZGNR (assuming collapsed band 
gap) at the Dirac point Ep = lO^^t, which is equivalent 
to the presence of disorder at the interface. 

Furthermore, the mode mixing is stronger, with larger 
off-diagonal t-matrix elements, for "lattice-unmatched" 



a, ^ 
g 2 



(a) 



• *^ "D 
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FIG. 13: (Color online) (a) The /- F characteristics, (b) trans- 
mission function, and (c) electric potential profile for 8-ZGNR 
(width ~ 1.5 nm; length « 4.4 nm) two-terminal device de- 
picted in Fig. 1121 The hopping parameter on both square and 
honeycomb lattice is tsi = t = 2.7 eV. In panels (a) and (c), 
dashed red line denotes AF insulating state of ZGNR, while 
solid blue line corresponds to metallic state with destroyed 
spin polarization and collapsed band gap. Panel (b) shows 
the transmission function in magnetic (red) and nonmagnetic 
(blue) states of ZGNR marked by points B and C in panel 
(a), respectively. Thin black line illustrates energy window 
f{E - eVds/2) - f{E + eVds/2) over which T{E, Vds) is inte- 
grated to get the current at the corresponding points B and 
C. The potential profile in panel (c) is plotted along the line 
drawn in the middle of ZGNR in Fig. [I2ja) . Thin vertical 
lines in panel (c) mark the boundaries of the extended device 
ZGNR -I- portion-of-metallic-leads, shown in Fig. 1121 across 
which the self-consistent voltage drop is calculated. 



leads in Fig. [TT] Therefore, the corresponding linear re- 

/ \ -1 

h 



sponse resistance R 



E 



of the device in 



Fig. [T2] (assuming collapsed band gap) is i? w 2.7 h/e'^ 
(the resistance quantum h/e^ is 25.8 kfl). On the other 
hand, the corresponding device with the same 8-ZGNR 
channel length {N^ = 37 or « 4.4 nm) and "lattice- 
matched" leads in the right column of Fig. [11] has R ~ 
2.29 h/e'^. The decay of the resistance with the 8-ZGNR 
channel length in the absence of defects, edge scatter- 
ing, impurities, and acoustic phonons (all of which were 
taken into account to extract the mean free path and 
contact resistance from R vs. plot in Ref. [l2 ) is due 
to reduced overlap of evanescent modes injected by two 
metallic electrodes.— For realistic contacts between vari- 
ous metals and graphene, one would also have to control 
the alignment of differing energy levels at the interface in 
order to reduce the effect of the Schottky barrier on the 
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contact resistancc i ' 

We choose the setup with "lattice-unmatched" leads 
in Fig. [12] as the device with lower contact transparency, 
to investigate their effect on the observability of voltage- 
driven collapse of edge magnetic ordering and the band 
gap corresponding to it. The current jump in Fig. [13] 
is much less pronounced than in Figs. [6] and [9] but it 
is still observable. This demonstrates that the value of 
the contact resistance itself, which is generated by dif- 
ferent transport mechanisms for devices in Figs. [7] and 
1121 is not the only reason for reducing the current jump 
at the threshold voltage. Moreover, despite poor con- 
tact between metallic leads and ZGNR, we still find hys- 
teretic behavior in Fig. [TS] which (assuming that it is 
not an artifact of self-consistent convergence algorithm 
for nonequilibrium electron densities, see Appendix |A| 
could also be used to confirm the band gap collapse in 
the nonequilibrium state of ZGNR. 

One can also compare the charge density redistribu- 
tion in Fig. [T2] and electric potential profile in Fig. [T3] 
to those in Sec. IV Bl for a two-terminal ZGNR device 
with different type of contacts. To ensure that elec- 
tric potential approaches the constant values in the bulk 
of the electrodes, a portion of the square lattice leads 
are attached to ZGNR to form an "extended device"— 
shown in Fig. [12] for which self-consistent calculations 
are performed. Thus, the charge transfer and poten- 
tial disturbance caused by ZGNR are screened off out- 
side the extended device region [voltage drop within the 
extended device region is enclosed by two vertical lines 
in Fig. [T^c)]. so that potential at the edges of the ex- 
tended device region matches to constant potential along 
the ideal semi-infinite electrodes. 



VI. CONCLUSIONS 

In summary, we predict voltage-driven nonequilibrium 
phase transition between magnetically ordered (Slater) 
insulating state and nonmagnetic metallic state in finite- 
length zigzag graphene nanoribbons. The ZGNR is at- 
tached to two metallic electrodes, where finite bias volt- 
age brings such two-terminal device into a nonequilib- 
rium steady-state with current flowing through it. The 
high density of states at the Fermi level, due to spe- 
cial topology of zigzag edges, results in instability when 
Goulomb interactions are taken into account, which is 
resolved through spin polarization around the edges 
in equilibrium. The spin-polarized state survives in 
nonequilibrium when the current flowing through ZGNR 
is small. However, at finite threshold voltage, the edge 
magnetic ordering is destroyed together with the band 
gap determined by the staggered potential of the magne- 
tization density profile. 

The abrupt jump of the current at the threshold volt- 
age, as one of our principal predictions, has unique and 
experimentally observable features: (i) current can in- 
crease by an order of magnitude when the bias voltage 



is tuned across the threshold voltage; (m) a hysteretic 
behavior in the I-V curve can occur around the win- 
dow of threshold voltages; (iii) the value of the thresh- 
old voltage increases with increasing nanoribbon length; 
{iv) the magnitude of the current discontinuity is reduced 
for poorly transparent contacts or metallic electrodes 
whose transverse propagating modes do not match well 
to ZGNR modes (propagating or evanescent) that can 
carry current. 

While these unique features can be tested with devices 
amenable to presently nanofabrication technology, they 
have not been observed in recent experiment o^^'^'^ on sub- 
10-nm-wide graphene nanoribbon-based field-effect tran- 
sistor devices which have displayed a sizable band gap 
for bias voltages up to 1 V. Thus, we delineate two pre- 
requisites for observing the collapse of the band gap and 
underlying magnetic ordering through the measurement 
of /- V characteristics of ZGNR two-terminal devices: 

1. The contact resistance between ZGNR and metallic 
electrodes should be kept low. This does not im- 
ply that the contacts must be perfectly transpar- 
ent. Our simulation results indicate that the I-V 
curve signatures of the band-gap collapse should 
be observable even when the contact transparency 
constitutes 20% from the ideal value, which corre- 
sponds to experimentally measured contact resis- 
tance of ~ 60 kfi in GNRFETs with top-deposited 
Pd electrodesJ^ 

2. The threshold voltage and the corresponding 
threshold current required to collapse the band gap 
of sub-lO-nm wide ZGNRs increase with increas- 
ing of nanoribbon length. The threshold source- 
drain voltage for sub-lO-nm-wide ZGNR which is 
« 6.6 nm long is w 0.44 V, and for ZGNR which 
is twice as long « 13.2 nm the threshold voltage Vt 
increases to « 0.8 V. At the same time, the shortest 
ribbon for which the experimental /- V curve mea- 
surements were performed^ was « 110 nm long 
and the applied source-drain voltage did not ex- 
ceed 1 V. This suggests that the threshold crite- 
ria has not been met. The decrease of the ZGNR 
length down to 10-20 nm range could result in ex- 
perimental observation of predicted current-flow- 
induced transition between spin-polarized and non- 
polarized ZGNR states. 

Thus, fabricating proposed devices — short graphene 
nanoribbon with atomically ultrasniooth zigzag edges^^ 
sandwiched between as transparent metallic contacts as 
possible — can be used to detect the presence of unusual 
s-p magnetism of carbon atoms along zigzag edges in un- 
ambiguous fashion and with all- electrical setup. If the 
predicted band gap collapse can be induced in such de- 
vices, the /-Incurve measurements will also probe aspects 
of spin dynamics in ZGNR. For example, if the current is 
turned off after ZGNR has been transformed from spin- 
polarized semiconducting to non-polarized semimetallic 
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state, the spin-polarized state is expected to be restored 
with some time delay Toff. That can be explained as fol- 
lows: the energy gain associated with small increment 
of the spin-polarized density at the zigzag edges is pro- 
portional to the spin density already accumulated at this 
edge. Because the spin density in non-polarized state is 
zero, no first order driving force is present to transform 
the system from non-polarized to spin-polarized state. 
The delay Toff can be measured as the time needed for /- V 
curve to change its character from conductive to highly 
resistive state. The time delay Ton in the onset of the 
band-gap collapse after current is turned on is expected 
to be much smaller than Toff. However, if Ton can be 
measured, information on electron velocity in the edge 
states, spin ordering, and the dependence of spin order- 
ing on temperature could be deduced, in principle. 
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APPENDIX A: SELF-CONSISTENT 
ALGORITHM FOR NONEQUILIBRIUM 
ELECTRON DENSITY IN ZGNR 

We compute the nonequilibrium solution using the 
voltage step Ay — 0.01 V. For the magnetically ordered 
ZGNR, the solution for V^J uses the solution for 
as the starting point. Calculations for the spin-polarized 
case start with vj"' = and proceed to the point Vt be- 
yond which the convergence towards the self-consistent 
solution cannot be obtained. Even though the inability 



to obtain convergence does not constitute a proof of the 
solution non-existence, we assume that Vt reached in this 
fashion is the threshold voltage at which spin polarization 
is destroyed. The Newton-Raphson algorithm described 
in Sec. IIIII guarantees the convergence towards the self- 
consistent solution, provided that sufficiently small por- 
tion of Sriin obtained from Eq. ([M]) is used to augment 
riin of the previous iteration. This self-consistent solution 
is not guaranteed to correspond to the lowest energy and 
depends on the initial electron density. Our algorithm 
monitors the convergence and adjusts the step along Sui^ 
vector to ensure that the riin is changed in a such way 
that the difference between riin and nout decreases during 
each iteration. 

When the step a becomes too small, i.e., the input 
density correction a x Jriin corresponds to the maximum 
potential variation of less than 1 meV that is still not 
small enough to make the difference between riin and 
riout in a current iteration less than in a previous one, 
the calculation stops. Even if the self-consistent solution 
can be obtained by further decreasing the step a, it will 
be highly unstable with respect to the external pertur- 
bations of the order of 1 meV. In a usual experimental 
environment this is equivalent to non-existence of spin- 
polarized solution. 

The non-polarized solution is obtained for sufficiently 
high applied voltage, such as Vds = 1.0 V, to ensure 
the solution stability. Then, the solutions for the lower 
voltages are computed using the higher voltage solution 
as the starting points. The calculations stop when the 
same criteria for the threshold voltage as in the case of 
spin-polarized solution is met. This can lead to well- 
pronounced hysteretic behavior of /- Incurves, depending 
on the contacts, as illustrated by Figs. [H^a), [HJa), and 
[IHa). 
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